% m_eval_1945_50_myopic.m
clear all; 
close all;
format shortG;
% *************************************************************************
% ****************** LOAD FUNDAMENTALS AND DATA  **************************
% *************************************************************************
load parameters.mat;                    
load calib_productivity_amenity.mat;   
load calib_setup.mat;           
load commuting_cost.mat;  
load calib_auxiliary.mat;        

load fundamentals_1950.mat;
load noagglomeration.mat; 

Rdata45 = Rdata_all(:,1); Rdata50 = Rdata_all(:,2); 
Ldata45 = Ldata_all(:,1); Ldata50 = Ldata_all(:,2); 
N = size(Ldata_all,1); 

% *************************************************************************
% ******************* PARAMETERS ******************************************
% *************************************************************************
alpha=alphaest;   
beta=betaest; 
theta50 = 0.90;
theta55 = thetavec(1); 

%% *************************************************************************
% *** OPTION VALUES OF AMENITIES AND PRODUCTIVITY EVALUATED USING 1945 ***
% *************************************************************************
XX50 = b50.*((Rdata45./Area).^beta); 
YY50 = a50.*((Ldata45./Area).^alpha); 
% ********************************************************************
% ******** SIMULATE FORWARD LOOKING MODEL ****************************
% ********************************************************************
R50_p = Rdata50; 
L50_p = Ldata50; 
Q50_p = Qdata(:,1); 

maxit=1e+6; 
ii=0; 
update=0.05;

while ii < maxit
    Ldif_p = L50_p-(1-theta50).*Ldata45;
    Rdif_p = R50_p-(1-theta50).*Rdata45; 
   
    w_p=(YY50).^(1/gammaL).*(Q50_p.^(-gammaH/gammaL)); 

    XX_temp = K50.*((repmat(Q50_p',N,1).^(-mu)).*repmat(XX50',N,1)).^(rho50/sigma50); 
    XX_temp = XX_temp./repmat(sum(XX_temp,2),1,N); 
    YY_temp = K50.*((repmat(Q50_p,1,N).^(-gammaH/gammaL)).*(repmat(YY50,1,N).^(1/gammaL))).^(rho50/sigma50); 
    YY_temp = YY_temp./repmat(sum(YY_temp,1),N,1);
    
    Lcom50_p = (1-theta50).*Lcom45+XX_temp.*repmat(Ldif_p,1,N);
    R50_n = (1-theta50).*Rdata45+sum(XX_temp.*repmat(Ldif_p,1,N),1)'; 
    L50_n = (1-theta50).*Ldata45+sum(YY_temp.*repmat(Rdif_p',N,1),2);
    
    H50_p = (1-psi).*hsupply45+(Q50_p.^eta).*Area; 
    Q50_n = (mu.*(w_p'*Lcom50_p)'+(gammaH/gammaL).*w_p.*L50_n)./H50_p; 
    Q50_n = (geomean(Qdata(:,1))/geomean(Q50_n)).*Q50_n; 
    
    R50_p_r = round(R50_p*1e+6)./1e+6; R50_n_r = round(R50_n*1e+6)./1e+6; 
    L50_p_r = round(L50_p*1e+6)./1e+6; L50_n_r = round(L50_n*1e+6)./1e+6; 
    Q50_p_r = round(Q50_p*1e+6)./1e+6; Q50_n_r = round(Q50_n*1e+6)./1e+6; 

    if min(R50_p_r == R50_n_r)==1 && min(L50_p_r == L50_n_r)==1 && min(Q50_p_r == Q50_n_r)==1
        break; 
    else
        ii=ii+1; 
        R50_p = (1-update).*R50_p+update.*R50_n; 
        L50_p = (1-update).*L50_p+update.*L50_n;
        Q50_p = (1-update).*Q50_p+update.*Q50_n; 
        if rem(ii,100)==1
        disp(['Largest Error is:',num2str(ii), ' ',num2str(max(abs(R50_p_r-R50_n_r))), ' ',num2str(max(abs(L50_p_r-L50_n_r))),' ',num2str(max(abs(Q50_p_r-Q50_n_r)))]);
        end
    end
end 

% *************************************************************************
% ************************** SAVE RESULTS *********************************
% *************************************************************************
result=array2table([ID, Rdata45, Rdata50, Ldata45,Ldata50, ...
    R50_p, L50_p, cbddist, Area],'VariableNames',{'geocode','pop45','pop50',...
    'emp45','emp50','pop50_memory','emp50_memory','cbddist', 'area'}); 
writetable(result,'../../output/quant/output_evaluations_1945_50_myopic.csv');
